#include <iostream>
#include <vector>
#include <cmath>        // std::abs
#include "CImg.h" 

using namespace cimg_library; 
using namespace std;

// @see https://www.omnicalculator.com/physics/harmonic-wave-equation#what-is-a-harmonic-wave
// @see https://www.tf.uni-kiel.de/matwis/amat/mw2_ge/kap_3/basics/t3_4_1.html
// @attention international SI units are used: meter, seconds etc.
/** 
    To determine the displacement, you can apply the following harmonic wave formula:

    y = A * sin[(2π / λ) * (x - vt) + Φ]

    where:

        y is the displacement of a given point along the wave,
        x is the position of that point (its distance from the source),
        t is the time point,
        v is the wave velocity,
        λ is the wavelength,
        A is the amplitude, and
        Φ is the initial phase of the wave.
 */ 

////////////////////////////////////////////////////////////////////////////////
typedef int position_t;  // data type of x and y coordinates along a surface
typedef double wave_t;  // data type of wave length and amplitude
////////////////////////////////////////////////////////////////////////////////
//~ // x- and y-value: can be used for row and column in jpg picture 
//~ // z-value/amplitude: can be used for rgb color for generating pictures
//~ struct Coordinate
//~ {
    //~ Coordinate( position_t xCord, position_t yCord, wave_t zCord=0 )
        //~ : x(xCord), y(yCord), z(zCord)
        //~ {;}
    
    //~ position_t x;
    //~ position_t y;
    //~ wave_t z;   // stores the amplitude of wave at a certain X/Y position
//~ };

////////////////////////////////////////////////////////////////////////////////
struct WaveSource
{
    WaveSource( position_t x_pos, position_t y_pos, 
                wave_t wLength, wave_t wAmplitude, 
                double wVelocity=0, double initPhase=0 )
        : xPos(x_pos), yPos(y_pos), 
          waveLength(wLength), waveAmplitude(wAmplitude), 
          waveVelocity(wVelocity), initialPhase(initPhase) 
        {
        }

    position_t xPos;      // x coordinate of the wave source
    position_t yPos;      // y coordinate of the wave source
    wave_t waveLength;    // the wave length
    wave_t waveAmplitude; // the amplitude of the wave
    wave_t waveVelocity;  // the velocity of the wave
    wave_t initialPhase;  // the initial phase of the wave in radians
    
    //unsigned int frequency;
};

////////////////////////////////////////////////////////////////////////////////
// consists of WaveSources and calculates wave amplitudes at certain positions
// @TODO: should we store the surface of waves in this class? TODO ?
class WavySurface
{
public:

    WavySurface( const position_t surfaceLength, 
                 const position_t surfaceWidth )
        : surfaceLength_(surfaceLength), surfaceWidth_(surfaceWidth)
    {
    }
    
    WavySurface( const position_t surfaceLength, 
                 const position_t surfaceWidth, 
                 WaveSource waveSource ) 
        : surfaceLength_(surfaceLength), surfaceWidth_(surfaceWidth)
    {
        waveSources_.push_back( waveSource );    
    }
    
    // adds a wave creating source to the surface
    void addWaveSource( WaveSource waveSource )
    {
        waveSources_.push_back( waveSource );
    }
    
    /** 
        To determine the displacement, you can apply the following harmonic wave formula:

        y = A * sin[(2π / λ) * (x - vt) + Φ]

        where:

            y is the displacement of a given point along the wave,
            x is the position of that point (its distance from the source),
            t is the time point,
            v is the wave velocity,
            λ is the wavelength,
            A is the amplitude, and
            Φ is the initial phase of the wave.
             
        x = 0.01 meter
        t = 1 second
        v = 0.005 m/s
        λ = 0.06024 meter
        A = 0.014 meter
        Φ = 0
            
        double y = 0.014 * sin( ((2*M_PI) / 0.06024) * (0.01 - 0.005*1) + 0);
        result: 0.00697469
     */
    wave_t calculateDisplacementAt(position_t x, position_t y, double timePassed)
    {
        wave_t displacement=0;
        for( vector<WaveSource>::const_iterator iter = waveSources_.begin();
             iter != waveSources_.end(); ++iter )
        {
            // get distance from wave source
            double distanceFromSource=getDistance(x, y, (*iter));    
            //cerr << "distance from source: " << distanceFromSource << endl;        
            displacement += ( *iter).waveAmplitude * sin( ((2*M_PI) / (*iter).waveLength) * (distanceFromSource - ((*iter).waveVelocity * timePassed)) + (*iter).initialPhase );
        }
        return displacement;
    }
    
private:
    
    inline double getDistance( position_t x, position_t y, const WaveSource& wSource )
    {
        // use pythagoras to calculate distance
        position_t xLen= abs(x - wSource.xPos);
        position_t yLen= abs(y - wSource.yPos);
        //cerr << "xLen: " << xLen << "\nyLen: " << yLen << endl;
        return sqrt( pow( xLen, 2 ) + pow( yLen, 2 ) );
    }
    
    const position_t surfaceLength_; 
    const position_t surfaceWidth_;
    vector<WaveSource> waveSources_;
};

// typedef Coordinate Pixel; 

////////////////////////////////////////////////////////////////////////////////
int main(void)
{
    // length and width of wavy surface (=picture display)
    //~ position_t dspWidth=1024;
    //~ position_t dspHeight=768;
    position_t dspWidth=1920;
    position_t dspHeight=1080;
    
    //create a wave generating source at position (0,0)
    //wave_t waveLen=20; // how long is one wave from crest to crest?    
    wave_t amplitude=30; // how high/low is wave?
    //wave_t velocity=0;  // how fast travels wave?
    WaveSource waveSource00(400, 100, 60, amplitude );
    WaveSource waveSource01(dspWidth/2, 3*dspHeight/4, 100, amplitude);
    WaveSource waveSource02(dspWidth, dspHeight, 40, amplitude);
    WaveSource waveSource03(-dspWidth, -dspHeight, 40, amplitude);
    
    // create a wavy surface
    WavySurface wavySurface( dspWidth, dspHeight );
    wavySurface.addWaveSource( waveSource00 );
    wavySurface.addWaveSource( waveSource01 );
    //wavySurface.addWaveSource( waveSource02 );
    //wavySurface.addWaveSource( waveSource03 );

    double secondsPassed=0;
    CImg<unsigned char> image(dspWidth, dspHeight, 1, 1, 128);

    for( int itX=0; itX<dspWidth; ++itX)
    {
        for( int itY=0; itY<dspHeight; ++itY)
        {
          //~ wave_t displ = trunc( wavySurface.calculateDisplacementAt(itX, itY, secondsPassed)+2*amplitude );
          //~ cerr << "displacement: " << displ << endl;
          //~ image(itX,itY,0) = displ;
          image(itX,itY,0) = trunc( wavySurface.calculateDisplacementAt(itX, itY, secondsPassed)+2*amplitude );
        }                
    }
    
    // save image
    image.save("WavePattern.png");

    CImgDisplay display(image, "Waves");
    while (!display.is_closed())
    {
        display.wait();       
        image.display(display);
    }
    return 0;
}
